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charge correction factors, the bandgap problem and the ab initio dielectric constant. 
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The errors arising in ab initio density functional theory studies of semiconductor point defects using the 
supercell approximation are analyzed. It is demonstrated that a) the leading finite size errors are inverse linear 
and inverse cubic in the supercell size, and b) finite size scaling over a series of supercells gives reliable isolated 
charged defect formation energies to around ±0.05 eV. The scaled results are used to test three correction 
methods. The Makov-Payne method is insufficient, but combined with the scaling parameters yields an ab initio 
dielectric constant of 1 1.6±4. 1 for InP. T point corrections for defect level dispersion are completely incorrect, 
even for shallow levels, but re-aligning the total potential in real-space between defect and bulk cells actually 
corrects the electrostatic defect-defect interaction errors as well. Isolated defect energies to ±0.1 eV are then 
obtained using a 64 atom supercell, though this does not improve for larger cells. Finally, finite size scaling of 
known dopant levels shows how to treat the band gap problem: in <200 atom supercells with no corrections, 
continuing to consider levels into the theoretical conduction band (extended gap) comes closest to experiment. 
However, for larger cells or when supercell approximation errors are removed, a scissors scheme stretching the 
theoretical band gap onto the experimental one is in fact correct. 

PACS numbers: 61.72.Bb 71.15.Dx 71.55.Eq 61.72.Ji 



I. INTRODUCTION. 

Understanding the properties of point defects and dopants is 
of key importance in studying the electrical and optical prop- 
erties of semiconductors. While various experimental tech- 
niques have been developed over the last half century it is 
only in recent years that they have started to be matched by 
accurate first principles computational techniques. Develop- 
ments in computing power have now made ab initio Density 
Functional Theory 1 (DFT) one of the most versatile atomic 
scale tools available for the investigation of defect properties 
in semiconductors and insulators. The key quantity to calcu- 
late is the defect formation energy 

E% = Ej {defect 9 ) - Ej(no defect) + ju,-n, - q (e v ± e F ) 

(1) 

where ^(defect) and Ej (no defect) are the total energy of 
the supercell "c" with and without the defect, (or charge q,) 
calculated using the same values of planewave cutoff, k-point 
grid, etc, to make use of the cancellation of errors. The defect 
is formed by adding/removing n,- atoms of chemical poten- 
tial fLj. £f is the Fermi level, measured from £,,, the valence 
band edge (VBE). Almost all properties of a defect can be 
derived from variations in and differences between formation 
energies. The method is very powerful, but critical limitations 
remain, two of the most important being the relatively small 
number of atoms which can be treated and the effect of the 
approximations, such as the Local Density and Generalized 
Gradient Approximations (LDA and GGA) required to solve 
the DFT itself. These treat quantum many body correlation 
and exchange effects incompletely, which in the case of semi- 
conductors and insulators results in a roughly 50% underesti- 



mation of the bandgap. This in turn has severe consequences 
for the calculation of defect transfer levels, the values of Ep at 
which the most stable charge state of the defect changes, given 
by the difference in E^ between the two states. The positions 
of the transfer levels govern whether the defect will be a single 
or multiple acceptor or donor, with levels deep inside the gap, 
or with shallow levels near to the band edges. Since the pre- 
dicted band gap differs so strongly from the experimental one 
it is very hard to map the calculated transfer levels onto the 
experimental gap and hence predict the electrical properties 
of the material. 

Meanwhile, the small number of atoms involved (100s or 
1000s) means that the boundary conditions become very im- 
portant. One of the most common approaches is to use Peri- 
odic Boundary Conditions (PBCs) together with a plane wave 
basis set 2 -. A supercell containing the defect in question is 
repeated periodically throughout space. The cell boundary 
thus looks bulk-like, rather than being a vacuum as with open 
boundary conditions. However, it also means that the defect 
interacts with an infinite array of images of itself seen in the 
PBCs. This alters E^ , making it (and most other defect prop- 
erties) supercell size dependent. The "true" defect properties 
are only recovered in the limit of an infinitely large supercell, 
equivalent to the limit of an isolated defect. This problem 
is particularly severe in the case of charged defects, where the 
Madelung energy becomes infinite if the charge is not neutral- 
ized using a uniform jellium background 3 . Even with jellium, 
the calculated formation energies can be wrong by several eV 
in supercells of the order of 10s or 100s of atoms, and we have 
previously shown 4,5 that finite size errors on this scale can 
even arise for neutral defects. Various authors, 6,7 ^ have at- 
tempted to create corrections schemes to estimate and remove 
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these errors, the most widely known being that of Makov and 
Payne 8 . Although these corrections are often used their accu- 
racy has been strongly questioned, with several studies sug- 
gesting that they are not reliably enough for regular use 6 - 910 . 

We previously 4 - 5 suggested that the supercell size errors can 
instead be eliminated by calculating the same defect proper- 
ties in a series of supercells of different sizes but the same 
symmetry and then finite size scaling the results to recover 
those of the infinite supercell. We found that varies with 
the supercell size, L, as 

£ d c (L) = £ d °° + a i L 1 + a n L-" (2) 

where a\, a„ and are fitting parameters, E°f being the 
finite size scaled formation energy corresponding to an in- 
finitely large supercell. The linear term has been discussed 
many times previously, first by Leslie and Gillian 3 . For neu- 
tral defects we found the correct value for n to be 3. This is 
actually very intuitive: most sources of error should vary with 
either the supercell size L (the defect-defect image distance) 
or with the cell volume L 3 (proportional to the jellium charge 
density, the number of atoms, the number of electrons, etc). 
Terms scaling with the surface area, 6L 2 , seem unlikely to be 
dominant. 

Here, two further sources of error must also be consid- 
ered. Firstly, since the electrostatic potential in a supercell 
with PBCs is only defined up to a constant, the zero on the 
energy scale must be chosen arbitrarily in each calculation. In 
the case of most pseudopotential codes (including the one we 
use) this occurs as an implicit average over values appropriate 
to each atom species in the supercell, weighted by the num- 
ber of atoms of each species. This means that the numerical 
value of Ef entering Eq [2 changes with the contents of the 
cell, leading to an additional finite size error. If the number of 
defects per supercell is constant then this error decreases with 
the number of atoms in the cell - essentially with the volume of 
the cell, L 3 . Hence this error is completely taken care of in the 
infinite supercell limit of our finite size scaling scheme. For 
individual supercells, Van de Walle and Neugebauer 10 suggest 
correcting the error by re-aligning the potential in the defect 
cell to that of the bulk, using its real-space value at some cho- 
sen point in a bulk-like region far from the defect. (We here 
use the point furthest from the defect in the unrelaxed cell.) 

Secondly, additional errors come from the dispersion of the 
defect levels introduced by overlap between the defect state 
wavefunctions and their PBC images. It has been suggested 1 1 
that this artificially raises E^ when k-points other than just 
the r point are used. It is suggested 11 that E% should then 
be shifted by q x (g£ — e|f ) , where e£ and £p are the val- 
ues of the Kohn-Sham level corresponding to the defect state 
calculated in the defect cell at the Y point and averaged over 
the sampled k-points respectively. The assumption is that the 
value of the defect level is correct at the Y point, so the differ- 
ence between that and the k-point averaged value should be 
removed. It has been shown by Hoglund et al. 12 that this is 
completely incorrect for the example of the phosphorus anti- 
site in GaP. By plotting the "bandstructure" of the defect level 
in different sized supercells it was shown that the defect level 
in the smaller cells is more-or-less correct when averaged over 



the sampled k-points, but much too low at the Y point. The 
same is also true for the As vacancy on the GaAs(llO) sur- 
face, for example 13 . Van de Walle and Neugebauer 10 instead 
point out 10 that in this respect there is a fundamental differ- 
ence between deep levels such as these and shallow defect 
levels. They suggest that the correction should only be ap- 
plied when evaluating transfer levels for shallow donors and 
acceptors. 

In the current paper we will show in section [TO] that Eq [2] 
with n — 3 also holds for charged defects, so that finite size 
scaling can be used to produce fully finite-size corrected de- 
fect formation and other energies, with well defined error bars 
and uncertainty. To do this we will study 1 1 example defects 
in the zinc -blend structured III-V semiconductor InP. These 
are chosen to include all types of native defects (vacancies, 
antisites and interstitials) as well as some common dopants 
at both substitutional and interstitial sites. Each is studied in 
one charge state only, usually the one that previous studies*^ 
suggest it has over the majority of the bandgap. The specific 
choices have been made to include all non-zero values from 
-3 to +3. 

These results will also enable us in section HV1 to perform 
the most objective and comprehensive reliability test we are 
currently aware of on other, computationally cheaper, correc- 
tion schemes. (Previous tests rely on only one or two - usually 
rather simple - examples, and do not generally have reliable 
isolated formation energies to compare with.) We will test the 
Makov-Payne scheme, potential re-alignment and dispersion 
corrections for shallow levels. In section ITVDI we will also 
derive an ab initio dielectric constant for InP by combining 
the scaling results with those of the Makov-Payne correction 
scheme. Finally, in section[V]we will use finite size scaling to 
provide the first clear-cut answer to the problem of mapping 
LDA or GGA transfer levels onto the experimental bandgap. 
Computational details are in the next section, and in section 
IVII we will conclude. 



II. COMPUTATIONAL DETAILS. 

We use planewave ab initio DFT 1 within the lo- 
cal density approximation (LDA) together with ultrasoft 
pseudopotentials 15 (US-PP) using the VASPcode 16 . Since we 
expect (at least) a 3 parameter fit we need at least 4 supercells. 
These must all be of the same symmetry since the errors scale 
differently for different symmetries. We choose simple cubic 
supercells containing 8, 64, 216 and 512 atoms. It would be 
preferable to replace the 8 atom cell by the 1000 atom one, but 
our computing resources are currently insufficient for k-point 
converged calculations with 1000 atoms. On the other hand, 
we previously found that, somewhat surprisingly, the 8 atom 
supercell is good enough in most cases: formation energies in 
this cell usually lie very close to the scaling curves, providing 
satisfactorily small error bars on the scaled values. Similarly, 
memory limitations force us to treat the indium 4d electrons 
as core, even though they are comparatively shallow: about 
14.5 eV below the VBE. This leads 5 to errors of up to ~ 0.5 
eV, but these are essentially supercell size independent. They 
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can easily be estimated in, say, the 64 atom cell and added 
back onto the scaled at the end. Our optimized LDA lat- 
tice constant using these chosen pseudopotentials 17 is 5.827 
A and the band gap is 0.667 eV, compared to 5.869 A and 
1.344 eV in experiment. We use jip = 3.485 eV and jii n = 
6.243 eV, corresponding to stoiciometric conditions, together 
with jU z „ = 1.891 eV, jU S i = 5.977 eV and jU s = 4.600 eV. For 
the 64 atom cell a planewave cutoff energy of 200 eV and a 
Monkhorst-Pack 4x4x4 k-point grid 18 was previously found 17 
sufficient to restrict errors to 0(0.01 eV) or less. When analyz- 
ing the errors arising from the supercell approximation itself, 
non-finite size dependent errors 5 (from the In pseudopotential, 
planewave cutoff etc) are not a problem. However, we do need 
to keep the k-point sampling errors down to at least the meV 
scale, since this convergence rate varies with supercell size. 
This is a much higher convergence criterion than is normally 
practical, necessary or even meaningful, and it is the reason 
that we pick only a limited number of example defects for this 
study. This convergence level can be achieved 5 by using the 
average 
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weighted by the number of points in the full Brillouin zone, 
where E^(N) is the formation energy calculated using an 
N x N x N Monkhorst-Pack k-point grid. The sum over TV 
is taken up to 12 in the 8 atom cell, 8 in the 64 atom cell and 
4, or for certain cases 6, in the 216 and 512 atom supercells in 
the unrelaxed geometries. (The weighted mean E^ converges 
much faster than the unweighted mean or the individual values 
Ej(N) themselves.) 

We present both non-relaxed (ions at ideal lattice sites) and 
relaxed calculations. No restrictions are placed upon the sym- 
metry of relaxations, but we do not allow atoms located on the 
surface of the cell to relax. The relaxation energy 
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(where E$ :Rx (N) and E$ M (N) are E$(N) with atoms at re- 
laxed and ideal positions respectively) converges faster with 
N than either E% :Rx (N) or E$' M (N). Hence we save compu- 
tational time by approximating the relaxed formation energies 
by 
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The relaxation energies used are weighted averages using 
6x6x6 and 8x8x8 k-point grids in the 8 atom cell, 2x2x2 
and (if the convergence is uncertain) 4x4x4 grids in the 64 
atom cell and 2x2x2 in the 216 and 512 atom supercells. For 
the latter cells we usually restrict the k-point grid to the ir- 
reducible Brillouin zone of the undisturbed bulk lattice. In 
other words, we use just the special k-point (0.25,0.25,0.25): 
the first Chadi-Cohen k-point 19 . This restriction is equivalent 
to assuming that the distortion in the bandstructure due to the 
presence of the defect is either localized (thus important only 
very near F) or symmetric. It introduces a small error whose 
significance again disappears in the large supercell limit. 



III. FINITE SIZE SCALING OF DEFECT FORMATION 
ENERGIES. 

A. Scaled formation energies for the example defects. 

Fig ^ shows the formation energy scaling for the 1 1 exam- 
ple defects in InP. The scaling curves using the uncorrected, 
as-calculated values are shown as solid lines in the figures 
(black in the online colour version). Their y-axis intersects 
give the Zs™ values listed in table [I] The curves also serve to 
predict the formation energy which would be expected in any 
finite sized supercell: for example the formation energies in 
the 8000 atom supercell are those at 1/L = 0.1 in the figures. 
We can estimate how accurate the E^ values are by adding the 
four dotted (black) curves shown for each example in Fig ^ 
in each of which one of the four data points has been omit- 
ted. (Note that for some cases the errors are so small that the 
dotted lines are hard to pick out, but they are still present in 
the figure.) The spread in y-axis intersects gives the error bars 
listed in the table. This is one of the particular advantages of 
using finite size scaling: it is possible not only to correct for 
the finite size errors themselves, but also to obtain a well de- 
fined uncertainty on the resulting energies - something other 
correction schemes can not provide. 

The errors obtained are on the 0.01-0.1 eV range or below, 
(smaller error are here rounded to 0.01 eV) and can doubtless 
be further improved if still large supercells are used. Note 
that, by construction, the errors which arise if only the 8, 64, 
and 216 atom supercells are used for the scaling are also on 
this 0.01-0.1 eV level (See tableHJ) 

The fact that such small error bars can be obtained indi- 
cates that a) scaling is a viable and practical approach to super- 
cell approximation errors, b) the k-point convergence is suf- 
ficient for our current purpose and c) our enforced use of the 
8 atom supercell is actually reasonable, for the same reasons 
described above and previously 5 . 



B. Form of the scaling. 

The choice of n — 3 again provides the best overall fit to 
the data, both for relaxed and non-relaxed calculations. Nor- 
malized % 2 tests 5 show that on average n — 2 provides fits 2.9 
times worse than n = 3 whilst n = 4 is 2.2 times worse. We 
note, however, that there are additional small (probably O(0. 1 ) 
eV or less) short-ranged errors present which decay exponen- 
tially with supercell size. These arise chiefly from the direct 
overlap of bound defect states with their PBC images and the 
resulting dispersion of the defect levels. In the case of relaxed 
energies some additional short ranged errors can appear be- 
cause defects in the 8 atom cell are only surrounded by 1 shell 
of relaxable atoms. The effect of this upon the form of the 
scaling can be seen in Fig [2] Here we show the scaling of the 
elastic contribution to the finite size errors. This is done by 
calculating formation energies in the 216 atom cell only, so 
that the electrostatic errors are essentially constant. The num- 
ber of shells of atoms permitted to relax around the defect is 
varied and the resulting formation energies are plotted against 



4 



T— tT— I 



r-H-H 




.1111 1 1 .. 


C): 


\ ~ 
•\ - 




: +3 

-In. 

i 

; ... i .... i . . 


■ 1 FT 



- -1 - 



I lllll I 1 1 -, II / Llllll I 1 r— |U 



\ \ : 
"•v \ - 




o _1 

Inverse supercell size (A ) (Top: number of atoms in cell.) 

FIG. 1: (Colour online.) Scaling of (x) unrelaxed and (+) relaxed formation energies. Curves are fits to Eq|5|with n=3. Solid (black) curves 
are fits to the four points as calculated (no corrections.) Dotted (black) lines each have one cell omitted for accuracy assessment. Scaling of the 
calculated values with various correction factors are shown for certain examples, as follows. Potential re-alignment: long dashed (red) lines 
in panels a) to g), and accuracy assessment for them: dotted (red) lines in panels a) to d). Dispersion corrections: short dashed (orange) lines 
in e) and f). Dispersion+potential corrections combined: dot-dashed (purple) lines in e) and f). First order (L _1 ) Makov-Payne corrections as 
dot-dot-dashed (blue) lines in panels g) to j). First+third order (L + L~ 3 ) Makov-Payne corrections: short dashed (pink) lines in h) to j). 
First order Makov-Payne+potential corrections combined: dash-dash-dot (green) lines in panel g). 



the inverse of the radius of the outermost relaxed atom shell. 
Hence the y-intersect corresponds to the formation energy ex- 
pected if an infinite number of shells are relaxed around the 
defect, but with the electrostatic errors inherent for the 216 
atom supercell. As expected, and as for the neutral defects 5 , 
the elastic errors are predominantly linear. Indeed, if the "one 
shell only" point from each curve is omitted then a linear fit 
works perfectly. (Solid lines in Fig [2]) The one shell only 
point corresponds to relaxations in the 8 atom cell, so we ex- 
pect that the elastic contribution to the supercell approxima- 
tion errors scales linearly with supercell size apart from some 
additional short range errors essentially only affecting the 8 
atom cell. 

These various short ranged errors have never-the-less only 
a very limited impact upon the final results, introducing only 
some additional scatter in the curves in Fig^ and hence lead- 
ing to larger error bars in some cases. They also lead to n = 2 
or n = 4 actually providing the best fit for some individual de- 



fects. However, in these latter cases the fitting with n = 3 is 
almost always a very close second. These problems can be 
overcome in a few years time once improved computing re- 
sources allow the study to be repeated using the 1000 atom 
supercell. For now we can still conclude that the elastic er- 
rors are essentially inverse-linear in supercell size, while the 
total formation energy errors (relaxed or unrelaxed) do indeed 
scale with the inverse-linear dimension and the inverse vol- 
ume of the supercell. 



IV. ASSESSMENT OF CORRECTION SCHEMES. 

In addition to the as-calculated formation energies, Fig ^ 
also shows the formation energy scaling using various cor- 
rection schemes. For clarity and space we do not show all 
possible corrections for all example defects, but results for 
all schemes are listed in tables HI and HTI All schemes recover 
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TABLE I: Scaled relaxed and unrelaxed (ideal lattice sites) formation energies with (ET ) and without potential corrections, for various 
example defects in InP. Note that the error bars are not actually symmetric: the widest has been listed in each case. £(L~ l ) and e(L~ 3 ) are 
ab initio values of the dielectric constant e for InP, calculated by comparing the Makov-Payne corrections of order L -1 and L~ 3 with the 
coefficients a\ and obtained from the scaling. e^(L _1 ) is the same thing calculated from the potential-corrected formation energies. All 
energies in eV, dielectric constants in units of the free space dielectric constant £q. 
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TABLE II: Assessment of correction schemes: Finite size errors (relative to the scaled values) are shown for the 64 atom supercell: S E is 
the error in the as-calculated formation energy, and 5 £+1+3 are the errors when Makov-Payne corrections are used to order L _1 and L -3 
respectively. S E ° A is the error when order L corrections are used, calculated with the ab initio dielectric constant evaluated from the results 
themselves (see text.) S E+k is the error when defect level dispersion is corrected for in the ionized states of the shallow donors and acceptors. 
Columns df etc are the same as S E etc but electrostatic potential realignments added. The averages are of the absolute error values |5 £ |. All 
energies in eV. 
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0.16 


-0.61 -0.50 


-0.07 


0.11 


0.06 


0.15 


0.07 


0.26 


-0.90 


-0.83 


Average 


0.64 


0.13 


0.19 


0.14 


0.07 


0.66 


0.46 0.49 


0.47 


0.29 


0.42 


0.18 


0.07 


0.75 


0.53 


0.64 






Average over both relaxed and unrelaxed structures: 


0.56 


0.21 


0.31 


0.16 


0.07 


0.71 


0.50 


0.57 



the correct formation energy in the infinite supercell limit, 
but not all produce improvements over the uncorrected for- 
mation energies for smaller supercells. This is shown in table 
im which lists the residual errors (relative to the infinite super- 
cell limit) when the corrections are applied in the 64 atom cell. 
The uncorrected 64 atom cell formation energies contain av- 
erage errors of about 0.5-0.6 eV, while using the potential re- 
alignment scheme produces errors of around 0. 1 eV. Makov- 
Payne does much worse (average errors around 0. 1 -0.4 e V, but 



often much larger) and the dispersion "corrections" produce 
errors which can be even larger than those in the uncorrected 
formation energies. 

A. Potential Realignment 

The potential re-alignment scheme is illustrated by the 
long-dashed (red) curves and points in Figs^a) to g). Even by 
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FIG. 2: (Colour online.) Scaling of the elastic contribution to the 
finite size errors in defect formation energies. Formation energies 
in the 216 atom shell are plotted versus the inverse of the radius of 
the outermost shell of atoms permitted to relax. ZnjJ (□), SiJ" n ' (y), 
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(red) lines: linear fits with the 1 shell only point omitted and included 
respectively. Dotted (blue) lines: quadratic fits to all points. 



the 64 atom supercell the values are very good indeed. How- 
ever, a lot of additional scatter is introduced into the corrected 
formation energies, . Indeed, the average errors relative 
to (the uncorrected) E°f do not shrink at all with increasing 
supercell size: 0.07 eV in the 64 atom cell, 0.10 eV in the 216 
atom cell and 0.09 eV in the 512 atom cell (see table[n}. This 
leads to wide error bars if the * values are scaled to give 

the infinite supercell limit, E°£'^ . We have derived scaling er- 
ror bars by the same technique described above. The resulting 
E7^ values and error bars are listed in Table [I] although the 
(red) dotted curves with data points are only shown in Figs^ 
a) to d). We find error bars of up to ±0.78 eV, average ±0.24 
eV. This means that potential realignment is a useful correc- 
tion for the results from individual supercells, but should not 
be used if more accurate results or defined error bars on re- 
sults are required. In that case non-realigned values should 
be scaled. These error bars are certainly too large to provide 
a basis for analysis of other correction schemes. The reason 
is that the correction scheme, good though it is, is not actu- 
ally complete or correct. Even in the largest supercells, the 
point furthest from the defect is not bulk-like, as the scheme 
assumes, resulting in either an over estimate or an under esti- 
mate, depending upon the specific conditions. 



B. Dispersion Corrections 

The dispersion correction scheme is illustrated for shallow 
donors and acceptors in Fig ^J) and f) and in Table Q] both 
with (short dashed, orange curves) and without (dot-dashed 
purple curves) potential alignment. Although the acceptor 
states fare better than the donors, the "corrected" values are 
always worse than those using only potential re-alignment, 
and usually worse than even the uncorrected formation en- 
ergies. Clearly, even for shallow defect levels, which follow 10 



closely the VBE or conduction band edge (CBE), still pro- 
duces worse formation energies than £q . 



C. Makov-Payne Corrections 

Fig[Dg) shows the first order L Makov-Payne corrections, 
with (dash-dash-dot, green) and without (dot-dot-dash, blue) 
potential re-alignment. When used together the two schemes 
usually produce a large over-estimate of the required correc- 
tion, (see columns 7 and 15 of table[ll]) almost as if using both 
corrections actually makes the same correction twice. Since 
the combination does so much worse than either technique 
alone there is no point going further with it. Instead, Figs^h) 
to j) show Makov-Payne corrections only, with formation en- 
ergies including both the order L 1 corrections (short dashed, 
magenta) and the order L _1 plus order L~ 3 corrections (dot- 
dot-dashed, blue). The order L _1 corrections work well in 
some cases (such as In^L when relaxations are omitted), but 
in most cases they are too large by a factor of about 1 A to 2, (as 
also noted by others 6 - 9 ) so that the "corrected" formation ener- 
gies are little better than the uncorrected ones. When the order 
L~ 3 corrections are added the correct formation energies are 
obtained in some cases, such as Vj 3 , but in other cases, such 
as Pj" n 2 , they help but are not sufficient. For other cases, such 
as Vp 1 , the corrections actually move the formation energies 
in the wrong direction. 

TablelTTIshows that the corrections are generally more likely 
to succeed for unrelaxed formation energies which is to be ex- 
pected since the electrostatic monopole terms are not the only 
ones to scale as L : the elastic errors do too. This means that 
even in principle the Makov-Payne corrections are only use- 
ful for non-relaxed formation energies, which are rarely the 
interesting ones. Besides this, the corrections also do better 
for more highly charged defects. This confirms that one of 
the problems is that they do not take into account the various 
other error terms which depend upon supercell size but not on 
charge state. These errors mostly have to do with the spuri- 
ous defect level dispersion introduced by the PBCs. Although 
the direct contributions of these are exponentially decaying 5 , 
their effects can still be seen in supercells on the scale of 10- 
100 atoms. Indeed the actual band width can be on the order 
of, for example, 0.5 eV and 2 eV in the 64 atom and 8 atom 
supercells 12 and remain significant even beyond that. Indirect 
dispersion effects can also be very important: for example, in 
a partially filled, erroneously dispersed defect level only the 
lower part will be filled, leading to too low a value for E^ . 
Worse happens if the defect level lies outside the band gap, 
either because it genuinely does or because the supercell is 
too small. This can lead to strong linear terms in the super- 
cell size errors even for neutral defects 5 : a neutral defect can 
behave as, say a -1 charged defect with (to a first approxima- 
tion) a +1 charged jellium background. This is not limited 
to neutral defects: a calculation for a defect anticipated to be 
in a +2 charge state (with a -2 charged jellium) could end up 
behaving more like a +3 charge defect with a -3 charged jel- 
lium. If the defect level moves outside the bandgap at certain 
k-points only it can lead to a linear error term which is not 
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even proportional to the square of an integer charge. Over- 
all, even the leading linear error term may be very different 
from that predicted by Makov and Payne's corrections. Un- 
fortunately, beyond noting that things get better on average 
for larger charge states and for non-relaxed calculations there 
seems to be no a priori method for determining whether the 
corrections will make things better or worse in a specific case. 
They are thus of little practical help, since they do not take into 
account enough of the specific behaviour of individual defects 
and materials. Indeed, it seems unlikely that any such highly 
generalized model for prediction of finite size error correction 
factors will ever fully succeed. 

D. Calculating the ab initio dielectric constant. 

Makov and Payne predicted that the two leading terms in 
the errors should be linear and cubic and our results show that 
they were correct in that respect. Their "corrected" formation 
energy takes the form 

E±mp(L) = E£ MP -h{EL)- l -h (eL) - 3 (6) 

where e is the dielectric constant, k\ = q 2 a/2 and k$ = 
2nqQ/3. (q is the charge of the defect, a is the Madelung 
constant for the supercell and Q is the quadrupole moment of 
the defect.) Comparing this with Eq|2]we find a„ = —k„/e. If 
we assume that the scheme is correct after all, then, q being 
known and Q having been calculated from the charge density, 
the only variable is the dielectric constant e. We can then 
use the correction scheme together with the scaling results to 
derive an ab initio value of £. This can be done twice for 
each defect, as shown in table|l] We find a wide scatter in the 
results, reflecting the wide variations in the effectiveness of 
the corrections. Indeed, the values of e obtained are com- 
pletely crazy when order L~ 3 corrections are used, as these 
are much more sensitive to short range effects and other er- 
rors. This, again reflects the fact that the situation described 
by Makov and Payne was highly idealized and ignores too 
many of the details of the charge distribution around specific 
defects. Never-the-less, the averaged values e from the or- 
der L 1 corrections are reasonably good. The most physically 
correct approach is to use the unrelaxed formation energies 
only, (with no elastic effects); indeed the values derived using 
relaxed values, third order corrections or Makov-Payne plus 
potential re-alignment make little sense (see|IJ. From the first 
order non-relaxed curves we obtain a dielectric constant of 
17.5±19.0. (The error bar is the standard deviation from the 
average.) Numerical problems with the Vp value have made it 
rather unreliable, and very different to the the others. Omitting 
it gives the perhaps more consistent value of 1 1 .6±4. 1 . These 
values compare to 9.6 in experiment or 10.7 calculated^ us- 
ing more traditional ab initio DFT-LDA techniques. 21 We thus 
obtain a fairly reasonable estimate of e as a free side-effect of 
performing accurate defect calculations - an interesting alter- 
native to the traditional calculations methods. The uncertainty 
in the value obtained is obviously rather large, but should im- 
prove if more defects in more charge states are included in the 
average. 



The order L Makov-Payne corrections do improve using 
this new value for the dielectric constant, see columns 5 and 
13 of table|J However, some individual values still have errors 
of up to 0.3-0.4 eV, and there is still no way to know when 
the corrections are making things better and when they are 
making things worse, so from a practical point of view the 
Makov-Payne scheme is still not reliable enough for accurate 
calculations. 



V. THE BANDGAP PROBLEM 

We now turn to the band gap problem and the issue of how 
to map calculated transfer levels onto the experimental gap. 
In practise several alternative - and essentially incompatible - 
methods are normally used. 

1) The Extended Gap Scheme: align the theoretical and ex- 
perimental VBEs and start plotting defect transfer levels from 
there, continuing past the theoretical CBE until one reaches 
the experimental one. In the section of the thus plotted "band 
gap" which lies above the theoretical CBE one automatically 
includes calculations in which supposedly localized, defect- 
bound electrons are in reality located in delocalized conduc- 
tion band states. The properties of the defect itself (transfer 
levels and local relaxed structure etc) re-enter primarily via 
hybridization of the conduction band states with the localized 
defect states, though this hybridization becomes smaller as the 
supercell size grows. 

2) The Scissors Scheme: align both the theoretical VBE 
and CBE with their experimental counterparts, performing a 
"scissors" operation to stretched out the theoretical gap states 
over the experimental gap. The manner in which this scissors 
operation should be done is not uniquely defined. A common 
option is to place acceptor levels the same distance above the 
experimental VBE that they appear above the theoretical VBE 
in calculations, and donor levels the same distance below the 
experimental CBE that they appear below the theoretical CBE 
in calculations. A better alternative is to actually examine the 
form and symmetry of the defect states themselves, and see 
whether they hybridize more strongly with host states near 
the CBE or with states near the VBE. If they hybridize most 
strongly with VBE states then they should be plotted the cal- 
culated distance above the (experimental) VBE, and if they 
hybridize most strongly with CBE states they should be placed 
the calculated distance below the CBE. 

3) The Reference Level Scheme. The basis of this scheme 
is rather different: the transfer level for the defect of interest 
is calculated, together with that of a similar reference defect 
for which the experimental value of the transfer level is well 
known, both done to the same level of accuracy. The differ- 
ence between the experimental and calculated levels of the 
known defect is subtracted from the calculated value of the 
new defect, so that the new level is only found relative to the 
old one. This idea is not without practical merit, but is very 
empirical. Its accuracy depends critically upon the choice of 
an appropriate reference defect, which must be as similar to 
the new one as possible, so it will not be discussed further 
here. However, it has an occasionally used ab initio variant, 
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FIG. 3: (Colour online.) Scaling of transfer levels for simple donors and acceptors calculated using: a) & b): neutral bulk as reference, with 
no corrections, c) & d): neutral bulk reference with potential corrections, e) & f) neutral bulk reference with dispersion corrections, g) & h) 
charged bulk as reference, with no corrections. Left panels: the dopant levels themselves. Right panels: the double donor or double acceptor 

Zn'^and 



levels, which should lie outside the bandgap. Using LDA: Sp" /0 and Sp" 2/+ (D, green), Si^fand Sif n 2/+ (0, red), Strand Sip /_2 (>, blue), 
Znj(~ 2 (A, pink). Using GGA: Sp'°(0) an d Zn^(v)- In g) & h) the smaller symbols with dashed lines show the acceptor-type levels relative 
to the experimental CBE rather than the LDA one. On all panels: the dotted lines are (in order of increasing energy) the VBE and CBE from: 
GGA (panel a)), LDA and experiment. The error bars shown have been constructed as described in section UTTI thoueh the dotted lines are 
omitted for clarity. 



which will be discussed: 

4) The Charged Bulk Reference Scheme: the reference state 
is not that of another defect, but is either the VBE or the CBE, 
meaning that a charged bulk total energy appears in Eq [2 
rather than a neutral one. In principle this provides an alterna- 
tive route around the band gap problem. (Details below.) 

Obviously, none of these schemes is fully correct, since the 
LDA/GGA bandgap problem is a fundamental one, but the 
important practical question of which approach comes closer 
to giving the correct physical picture remains unanswered. In 
principle it can be answered by examining various experimen- 
tally well known defect levels. The exact location of most 
native defect levels is rather hard to measure to a sufficiently 
high accuracy to answer this question, but many simple donor 
and acceptor levels are known very accurately. We will use 
the 0/- acceptor level of Zni n , which in experiment lies 0.035 
eV from the VBE, and the +/0 donor levels of Sp and Sii„, 
which experiment finds about 0.006 eV from the CBE. We 
will also add in the 0/- transfer level of Sip, which would be 
a simple acceptor if Sii„ had not been the more stable site for 
Si in InP. This gives us an example of a donor and an acceptor 



on each sublattice, so that all bonding and band hybridiza- 
tion possibilities are represented. Unfortunately, calculations 
of these levels in finite sized supercells in the 100-200 atom 
range have never produced a clear answer to the question, so 
we will use finite size scaling to correct for the supercell ap- 
proximation errors. The results are shown in Fig[3]a) using as 
calculated values, c) adding in potential corrections and e) us- 
ing dispersion corrections. (Van der Walle and Neugebauer 10 
suggested that dispersion corrections should still be correct for 
shallow transfer levels.) The results using as calculated trans- 
fer levels and potential corrected ones are very similar. The 
dispersion corrections, on the other hand, are clearly com- 
pletely incorrect: they place both acceptor and donor levels 
in the midgap for most practical supercell sizes, whether the 
potential corrections are added (not shown) or omitted (as 
here). Meanwhile, in Fig[3]b), d) and f) we also show the 
second donor/acceptor levels, Znj^ 2 , Sp 2/+ , etc., calculated us- 
ing the same correction schemes. Since these levels are never 
observed experimentally they must lie outside the band gap. 
Hence the VBE should lie between the double donor levels 
(right panels) and the single acceptor levels (left panels). Sim- 
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ilarly the CBE should lie between the single donor and double 
acceptor levels. In practise, these pairs of levels more or less 
coincide, doubtless a result of the remaining limitations in the 
use of DFT-LDA for semiconductor defects. Fortunately this 
still leaves us with a clear view of how to treat the band gap 
problem. 

In the 64 atom cell the donor (and double acceptor) lev- 
els lie roughly the experimental band gap (1.3 eV) above the 
VBE, while the acceptor (and double donor) levels lie on av- 
erage a little below the VBE. However, coming to the larger 
cells the donor levels fall and the acceptor levels rise. Finite 
size scaling places the acceptor levels Zn^and Sip"0.03 and 
0.01 eV above the VBE respectively, in rather good agreement 
with experiment. The single donor (and double acceptor) lev- 
els all scale to the theoretical CBE. 22 To be more specific, 
transfer levels calculated using LDA scale to the LDA band 
edges, while the Zn^~ and Sp /0 transfer levels calculated using 
the Perdew-Wang GGA 23 scale to the edges of the GGA band 
gap - Fig a). (The GGA CBE lies 0.2 eV below the LDA 
one when the lattice parameter has been optimized.) 

Hence, scheme 1, the extending gap scheme, is seen to be 
the most appropriate when only reporting uncorrected results 
from supercells of about 50-100 atoms. However, when the fi- 
nite size errors are removed (by scaling or by some other tech- 
nique) it becomes clear that the scissors scheme, scheme 2, is 
physically far more correct. Unsealed LDA or GGA results in 
supercells over a few 1000 atoms would also be best reported 
using the scissors scheme. For intermediate (100-1000 atom) 
supercells some kind of hybrid approach is required. The re- 
sult also indicates why the controversy has lasted so long: ul- 
timately the scissors scheme is correct, but this only shows up 
for very large supercells or with scalingi2i 

We now return to scheme 4), the charged bulk reference. 
This amounts to replacing the terms —Ej (no defect) and 
—q£F in Eq^by the term — Ej (no defect^), which is the to- 
tal energy of the bulk supercell c with — q extra electrons and 
neutralizing jellium. Fig |3] g) & h) show the transfer levels 
calculated like this, with no corrections terms. The donor lev- 
els behave in the same way as using Eq^in Fig[3]a), but the 
acceptor levels are less straight-forward. Using a charged bulk 
reference the levels come out relative to the CBE, rather than 
the VBE: they implicitly include the bandgap, which must be 
subtracted off again to place them on the same overall scale 
as the donor levels. This gives a "choice" for the value for 
the bandgap to subtract, which is how the potentially route 
around the band gap problem enters. Namely, if the 0/- trans- 
fer level emerges as, say, -0.5 eV, one could place it 0.5 eV 
below the experimental CBE, thus plotting the transfer levels 
over the experimental bandgap. (Small symbols and dashed 
scaling curves in Fig[3]g & h).) For the single acceptor lev- 
els this clearly does not work: although they land accidentally 
close to the VBE for smaller supercells they actually scale to 
the theoretical CBE, which is completely wrong. Instead, they 
should be placed below the theoretically CBE (large symbols, 
solid curves), where they scale to the VBE. Unfortunately the 
opposite is true for the double acceptor levels. These work out 
roughly OK if plotted relative to the experimental CBE - lying 
outside the theoretical band gap, even if still inside the exper- 



imental one - but using the theoretical CBE, (as required for 
the single acceptors) they lie inside the theoretical band gap, 
disagreeing with experiment. Hence using a charged bulk to- 
tal energy as the reference for charged defect calculations is 
not even internally consistent and the scheme is thus funda- 
mentally incorrect. 



VI. CONCLUSIONS 

In this paper we have shown that finite size errors in the 
supercell approximation scale with the linear dimension and 
with the volume of the supercell, and that finite size scal- 
ing the results from a series of supercells removes the su- 
percell approximation errors, leaving accurate information on 
isolated semiconductor defects, without the need for correc- 
tions. We also obtain error bars defining the uncertainly on 
the results obtained, and as far as we are aware this is the only 
method which is able to remove these errors in a controlled 
manner with defined uncertainty. We have demonstrated this 
using a variety of different types of defects with charge states 
ranging from -3 to +3 and find that it is possible to reduce for- 
mation energy errors from the 0. 1 -2 or so e V range of practical 
supercells down to the 0.01-0. 1 eV range or below - doubtless 
much lower if still larger supercells are used. By construction, 
errors on this scale also occur if only the 8, 64, and 216 atom 
supercells are used. 

We then used the scaled results for the first full reliability 
test of three correction schemes. We found that dispersion cor- 
rections are incorrect and Makov-Payne corrections are poor 
(with both the experimental and LDA dielectric constants,) 
though they did allow us to obtain a reasonable ab initio LDA 
dielectric constant of e = 1 1 .6±4. 1 for InP. On the other hand, 
the potential re-alignment scheme was found to be remarkably 
successful, removing much of the electrostatic defect-defect 
error as well, to leave average residual errors of about 0. 1 eV, 
from single calculations with supercells in the 64-512 atom 
range. 

This obviously raises the question of why the potential re- 
alignment scheme is so successful, when it does not set out to 
correct defect-image interaction errors at all! The fact that it 
produce similar (but more reliable) corrections to the Makov- 
Payne scheme suggests that it is some how dealing with the 
electrostatic errors anyway. We noted in section HV Al that the 
scheme assumes that the real-space potential at some point in 
the cell far from the defect is bulk-like, even though for practi- 
cal cell sizes it is not bulk-like at all. The resulting additional 
shift in this local real-space potential reflects the effects of 
the electrostatic defect-image interactions. Doing the poten- 
tial re-alignment in this way therefore fails to properly cor- 
rect the mismatch in the zeros of the energy scales between 
the bulk and defect cells, but the "error" in the re-alignment 
more or less corrects for the electrostatic errors arising from 
the PBCs. 

Finally, we have given the long awaited answer to the 
dilemma of how best to map LDA and GGA calculated defect 
transfer levels onto the experimental gap, and indicated why 
the issue was previously so hard to settle. The key result is 
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that the scissors method is physically more correct, though the 
extended gap scheme is best when reporting results from sin- 
gle supercells on the 1 -200 atom scale without finite size cor- 
rections. For uncorrected results from supercells over a few 
1000 atoms the scissors method is best, with a hybrid method 
needed in between. The best, of course, is to use the scissors 
scheme, with either scaled or corrected results, regardless of 
supercell size. The apparent success of the essentially incor- 
rect extended gap scheme for uncorrected results in manage- 
ably sized supercells is the basic reason for the debate lasting 
so long. 

This leads to another issue which is also apparent from our 
results: It is very dangerous to report calculations from sin- 
gle supercells without trying to estimate the errors contained. 
Quantitatively these can be ~l-2 eV or more, but we have 
cases here where conclusions are even qualitatively wrong in 
supercells up to and even including the 512 atom cell. For 
example 25 , comparing these results for P^ with those 5 for Pj^' 
we find that even at the CBE, the +2 charge state appears more 
stable than the +0 in all four supercells. The fact that it is ac- 
tually 0.19 eV less stable only emerges when the finite size er- 
rors are removed, either by scaling or (leaving residual errors 
from 0.05-0. 13 eV) by using potential realignment. Similarly, 
at the VBE, Vj„ and Inj^ are more stable than Vj" n ° and Injpj, 
respectively, in both the 64 and 216 atom supercells. The cor- 
rect stability order only appears in the 512 atom cell, (neu- 
trals more stable by 0.21 and 0.16 eV respectively,) and the 
correct order of magnitude for the difference (0.68 and 1.17 
eV) is only obtained by scaling. Another striking example is 
that, according to LDA in cells <512 atoms, p-type Zn doped 
InP - a material upon which much of current optoelectronics 
depends - should not be p-type at all! For the roughly stoicio- 
metric conditions of, say, Czochralski growth, LDA in the 64 
and 216 atom cells places Zn not as the shallow acceptor Zni n 
but as the interstitial Zn^p), where it is a deep double donor. 
Even in the 512 atom cell the two are degenerate, suggesting at 
best semi-insulating material. According to this Zn is only an 
acceptor for InP grown under strongly non-equilibrium condi- 
tions, such as with molecular beam epitaxy. However, Zn is a 
p-type dopant, even grown from the melt, and this fact can be 
predicted using LDA, but only for supercells of the order of 
1000s of atoms, or if the supercell size errors are removed - by 
scaling or otherwise. Doing this using potential realignment 
works for all these examples: even in the 64 atom supercell 
reasonable results can be obtained. However, caution should 
still be used: Firstly, for our examples it worked much better 
for the formation energies than for the shallow dopant transfer 



levels. Secondly, potential re-alignment makes Pj" n (correctly) 
more stable than P^ at the CBE in all but the 8 atom supercell, 
but if those corrected results are then scaled the wrong answer 
returns, with PJ" 2 more stable than Pj" n ° because of the large 
error bars found when scaling potential re-aligned energies. 

In short, it is essential, to estimate and report the finite size 
errors for each specific case when reporting supercell defect 
calculations. This is often omitted, or is only done using the 
unreliable Makov-Payne scheme. When it is done this is usu- 
ally by doing most calculations in a cell of, say, 50-200 atoms, 
and then repeating a few of them in a slightly larger cell. If 
the calculated results do not change much then they are con- 
sidered converged. However, even this should be done with 
extreme caution. Even with only a linear contribution, the fi- 
nite size errors in the 64 atom supercell are three times the 
difference between the 64 and 216 atom cell energies, the 216 
atom cell errors are still twice this estimate. Even the errors 
in the 512 atom cell are three times the difference between the 
216 and 512 atom energies. 

So, how should finite size errors within the supercell ap- 
proximation be treated? Ideally, using finite size scaling of 
otherwise uncorrected energies. This is, of course, costly in 
both human and computer time. The best alternative is sim- 
ply to use potential realignment in as large a supercell as time 
and resources permit. However, one should be aware that a) 
this does not help the elastic errors, b) potential realignment 
should not be combined with finite size scaling and c) there 
is no way to estimate the remaining errors or the reliability of 
the results. For our examples, the average errors using this 
method are ^0.10 eV, but with some examples up to 0.21 eV, 
and nothing to say that much larger errors will never occur. 
If the conclusions being drawn from a calculation are not ad- 
versely affected by uncontrolled errors of 0.1-0.2+ eV then 
this method is reasonably good. Otherwise, the only truly reli- 
able method of controlling the errors in the supercell approx- 
imation, and defining the uncertainly in the results, is finite 
size scaling. 
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tions or scissors shifts of the bandstructure, and so comes from a 
similar level of method to ours. 

The scaling of the single donor and double acceptor levels pro- 
duces rather wide error bars. The reason is that the particular 
Kohn-Sham level which is half filled in the charge state but 
empty in the +1 charge state lies above the CBE in all of these 
(finite) supercells, though it will re-appear inside the gap for large 
enough cells. As we showed previously 5 this tends to produce 
poor scaling in the charge state. 
J.P Perdew and Y. Wang Phys. Rev. B 13, 5188 (1976) 
The data shown in figure|3|were calculated using the In 4d elec- 
trons as core, but the conclusions are completely unaffected: cor- 
recting this error using In 4d valence calculations in the 64 atom 
cell, the curve for Zn^" moves down by about 0.003 eV, that for 
Sip'" moves up about 0.04 eV and those for Si^° and Sp" /0 move 
down about 0.08 eV. We also note here that the present results are 
for relaxed formation energies, but that we can reach exactly the 
same conclusions using non-relaxed formation energies. 
The numerical results reported in this paragraph include correc- 
tions for the use of the In 4d core pseudopotential. They are thus 
considered accurate to around 0.01-0.02 eV, not counting errors 
in the LDA itself. 



